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Abstract. A new approach to the problem of gravitational waves detection based on simultaneous 
timing of several pulsars and subsequent expansion of the post-fit timing data into components of 
different spectral kind (with different spectral indices) is proposed. Presence of a signal caused 
by stochastic gravitational waves in spectral components is tested with the two-point angular 
correlation function as proposed in the pulsar timing array. This new approach was applied to timing 
data of a few millisecond pulsars and allowed to detect a signature similar to one predicted for 
gravitational wave background at relatively high confidence level: correlation coefficient between 
experimental and theoretical two-point correlations p — 0.82 ± 0.07. 
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INTRODUCTION 

Fractional instability of some pulsars is comparable with the one of atomic standards 
and reaches the level 10^^^. Such a property allows to use pulsar timing for solution of 
various astronomical and metrological problems. In this paper a method of detection of 
the gravitational waves (GW) by pulsar timing array based on the two-point correlation 
function (eq. 1) is further developed [5, 6, 14]. 

m = \^^og{x)-'^+\+\d{xi (1) 

where x = ( 1 — cos 0)/2, 9 is angular distance between pulsars, 5 is delta-function. 

Many radio observatories started pulsar timing observational programs with the final 
aim of detection of the gravitational waves. By the year 2005 the Kalyazin Radio 
Astronomical Observatory (KRAO) of the Astro Space Center of the Lebedev Physical 
Institute already had 8-9 year time series of pulsar timing data of 6 millisecond pulsars 
(PSRs J0613-0200, J1640-I-2224, J1643-1224, J1713-I-0747, J1939-I-2134, J2145-0750). 
After publication of the paper [6] method of the two-point correlation function was 
applied to the data (pulsar post-fit residuals of time of arrivals (TOAs)) and showed 
that it can not be used directly since TOA residuals always display a presence of so 
called timing noise that prevents detection of the signature (1) clearly. 

However the additive noise presented as a phase variations in all frequency standards 
(atomic clocks or pulsars) has distinct features which allows to separate one kind of the 
noise from another. Each kind of the additive noise can be described in form /j„/(u", 
where n is spectral index of the power spectrum of the phase variations, (O is Fourier 
frequency, /?„ is the amplitude of the noise. Each spectral index n corresponds to noise 
caused by different physical phenomena [2]. E.g., n = corresponds to the measurement 
noise, n = 3 corresponds to the phase variations due to interstellar medium [2], n = 5 



corresponds to the variations caused by the primordial gravitational wave background 
(GWB) [1]. 

As an experience of frequency standards exploitation shows, noise with more steep 
spectrum begins to dominate at more long time interval since the amplitudes hn rapidly 
diminish with increasing of n. This property of the frequency standards can be used for 
separation of the additive noise into components with different spectral indices by an 
appropriate mathematical method. A possible presence of the GWB in the components 
can be detected with the two-point correlation function. The singular spectrum analysis 
(SSA) was used for expansion of pulsar timing data into components [4]. As fig. 5 from 
the paper [12] shows, the SSA-method indeed expands the input stochastic signal (TOA 
residuals) into components with different spectral indices. This property is particularly 
noticeable for low and high frequency components. 

OBSERVATIONS 

The pulsars PSR J0613-0200, J1640+2224, J1643-1224, J1713+0747, J1939+2134 and 
J2145-0750 were observed with the fully steerable RT-64 radio telescope of KRAO [8, 9, 
11]. The AS -600 instrumental facility of the Pushchino Radio Astronomy Observatory 
(Astro Space Center, Lebedev Physical Institute) was used for the registration [10]. The 
pulsar pulses were accumulated by a spectrum analyser in two circular polarizations, 80 
channels in each, with a channel frequency band 40 kHz. The observing sessions were 
conducted, on average, once every two weeks. The total signal integration time in each 
session was about 2 hours. 

The topocentric TOAs were determined by fitting the session-summed pulsar pulse 
profile into a reference template with a hight signal to noise ratio. We computed the 
barycentric TOAs, determined the TOA residuals and refined the pulsar timing parame- 
ters by minimizing the residuals by the least-squares method using the software Tempo 
[13]. The DD model [3] was used to refine and compute the pulsar orbital parameters. 
The astrometric, spin and orbital parameters of all six pulsars can be found in the paper 
[12]. Fig. 1 shows the post-fit TOA residuals of 6 pulsars. 

COIMPUTER SIMULATION 

For computer simulation the coordinates of above six pulsars were taken. A total of 
N — 50 generations were made in accordance with the number of TOAs. The amplitude 
of the GW was set to constant h=l, right ascension of the GW was uniformly distributed 
in the range 0^<a< 24^ and sine of declination of the GW was uniformly distributed 
in the range — 1 < sin 5 < 1 . Next, the angular coefficients were calculated on the basis 
of formula (1) from [5]. Fig. 2a shows results of the computer simulation. On the basis 
of the simulation one can conclude that with the M = 6 pulsars and N = 50 TOAs the 
two-point correlation function ^{9) is detected with the amplitude signal-to-noise ratio 
SNR ~ 10. 
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FIGURE 1. Post-fit timing residuals of six millisecond pulsars observed at the Kalyazin Radio Astron- 
omy Observatory. 
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FIGURE 2. a) Results of the numerical simulation of the two-point correlation function for M = 6 pul- 
sars and = 50 TOAs. Signal-to-noise ratio SNR ~ 10. b) Experimental plot of the two-point correlation 
function. Moving average by 2-5 points were applied. Correlation coefficient between theoretical curve 
(solid line) and experimental points p = 0.82 ± 0.07 



RESULTS 



For simplification of the computer modelling all pulsar data were binned and averaged 
in td = 40 days intervals. Gaps were filled by the linear interpolation of adjacent 
values. The common part of the observation interval (~ 6 yr) was taken only. Averaging 
were calculated at the different intervals IS'^ < A? < 60'^. The discussed effect reveals 



itself most clearly at /St ~ 40 — 50 . To make the experimental plot of C,{9) more 
dense in -domain the moving average by 2-5 points was applied to the plot. Fig. 2b 
shows experimental plot of C,{9). Correlation coefficient between theoretical curve and 
experimental points p = 0.82 ±0.07. For this correlation coefficient and the number 
of points M(M— l)/2 = 15, the probability that this correlation occurred by chance is 
Pr(p = 0.82 ±0.07, M = 15) < lO^^ [7]. 

CONCLUSION 

A new approach to the problem of GW detection with pulsar timing array is proposed. 
This approach is based on the fact that real physical time series always contain a 
superposition of noise signals of different kind (with different spectral indices). As a 
rule the signals have a characteristic amplitude and begin to dominate at characteristic 
intervals. This feature allows to expand the time series into components of different kind 
using an appropriate mathematical method. As such a method the singular spectrum 
analysis was used. Combination of two methods, SSA and pair angular correlation, 
allowed to detect in TOA residuals a signal (correlated amplitude ~ 0.5 mcs) which 
at statistically significant level produces similar with GWB effect. 
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